La presente analisi si propone di affrontare un problema di classificazione che ha come obiettivo discriminare tra individui fumatori e non sulla base di una serie di misurazioni mediche. Le tecniche di data mining applicate sono: alberi decisionali, KNN, regressione logistica e random forest.
I dati utilizzati possono essere scaricati al seguente link: https://www.kaggle.com/datasets/kukuroo3/body-signal-of-smoking
Il fumo è una delle più grandi minacce per la salute in quanto rappresenta uno dei maggiori fattori di rischio nello sviluppo di patologie neoplastiche, cardiovascolari e respiratorie. Secondo i dati diffusi dall’OMS è responsabile del decesso di almeno 6 milioni di persone ogni anno. Questo si traduce in ingenti spese sostenute dai servizi sanitari nazionali per le cure di pazienti affetti da patologie causate dal fumo.
L’esempio qui considerato è quello del National Health Insurance Service (NHIS), sistema di assicurazione sanitaria nazionale del Sud Corea. Nel 2014 il NHIS ha indetto una causa contro KT&G, il principale produttore di tabacco coreano, richiedendo un risarcimento di 53,3 miliardi di won, pari alla spesa totale sostenuta dal 2003 al 2013 per il trattamento di pazienti affetti da cancro ai polmoni e alla laringe, fortemente legati al fumo.
L’NHIS ha raccolto varie informazioni sui check-up svolti da un campione di cittadini assicurati dal servizio stesso. Queste sono costituite da dati identificativi dei pazienti, dai risultati della visita e dallo stato di fumatore.
La seguente analisi applicherà quindi diversi metodi di classificazione al fine di identificare lo stato di fumatore per pazienti di cui si conoscono solo valori biologici. Questo permetterebbe all’NHIS di avere una previsione su quanti tra i loro assicurati siano fumatori o meno.
Il dataset contiene 55692 osservazioni. La variabile target Smoking assume valore 0 nel caso di individuo non fumatore e 1 in caso di individuo fumatore. Il 63% degli individui è non fumatore, il 37% lo è. Sono inoltre presenti 26 variabili così definite:
smoke <- read.csv("smoking.csv")
smoke$gender<-as.factor(smoke$gender)
smoke$hearing.left.<-as.factor(smoke$hearing.left.)
smoke$hearing.right.<-as.factor(smoke$hearing.right.)
smoke$dental.caries<-as.factor(smoke$dental.caries)
smoke$tartar<-as.factor(smoke$tartar)
smoke$smoking<-as.factor(smoke$smoking)
Partendo dal dataset originale appare opportuno procedere subito all’eliminazione di alcune variabili, quali
ok <- ggplotly(ggplot(smoke, aes_string(x = "smoking", y = smoke$eyesight.right., col = "smoking", fill = "smoking"))+geom_boxplot(alpha=0.2)+ylab("eyesight.right"))
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ok2 <- ggplotly(ggplot(smoke, aes_string(x = "smoking", y = smoke$eyesight.left., col = "smoking", fill = "smoking"))+geom_boxplot(alpha=0.2)+ylab("eyesight.left")+labs(title="boxplot di eyesight.right sopra e di eyesight.left sotto"))
subplot(ok,ok2,nrows=2)
library(dplyr)
smoke<-smoke%>%dplyr::select(-c(ID,eyesight.left., eyesight.right., Urine.protein, oral)) #rimozioni variabili
str(smoke)
## 'data.frame': 55692 obs. of 22 variables:
## $ gender : Factor w/ 2 levels "F","M": 1 1 2 2 1 2 2 2 1 2 ...
## $ age : int 40 40 55 40 40 30 40 45 50 45 ...
## $ height.cm. : int 155 160 170 165 155 180 160 165 150 175 ...
## $ weight.kg. : int 60 60 60 70 60 75 60 90 60 75 ...
## $ waist.cm. : num 81.3 81 80 88 86 85 85.5 96 85 89 ...
## $ hearing.left. : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ hearing.right. : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ systolic : num 114 119 138 100 120 128 116 153 115 113 ...
## $ relaxation : num 73 70 86 60 74 76 82 96 74 64 ...
## $ fasting.blood.sugar: num 94 130 89 96 80 95 94 158 86 94 ...
## $ Cholesterol : num 215 192 242 322 184 217 226 222 210 198 ...
## $ triglyceride : num 82 115 182 254 74 199 68 269 66 147 ...
## $ HDL : num 73 42 55 45 62 48 55 34 48 43 ...
## $ LDL : num 126 127 151 226 107 129 157 134 149 126 ...
## $ hemoglobin : num 12.9 12.7 15.8 14.7 12.5 16.2 17 15 13.7 16 ...
## $ serum.creatinine : num 0.7 0.6 1 1 0.6 1.2 0.7 1.3 0.8 0.8 ...
## $ AST : num 18 22 21 19 16 18 21 38 31 26 ...
## $ ALT : num 19 19 16 26 14 27 27 71 31 24 ...
## $ Gtp : num 27 18 22 18 22 33 39 111 14 63 ...
## $ dental.caries : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 1 1 1 ...
## $ tartar : Factor w/ 2 levels "N","Y": 2 2 1 2 1 2 2 2 1 1 ...
## $ smoking : Factor w/ 2 levels "0","1": 1 1 2 1 1 1 2 1 1 1 ...
Innanzitutto si suddivida il dataset in training e test set. Il primo verrà usato per la fase di pre-processing e di scelta e valutazione dei modelli, il secondo per la classificazione delle nuove osservazioni e valutazione delle previsioni. Avendo una numerosità elevata di osservazioni si è scelto di effettuare una divisione 70% train - 30% test.
set.seed(123)
index <- sample(1:nrow(smoke),0.7*nrow(smoke))
train <- smoke[index,] #train: 38984 obs
test <- smoke[-index,] #test: 16708 obs
Dato che saranno valutate le performance di più modelli, si applica un’ulteriore divisione all’interno del set di training creando così un nuovo training set e un validation set.
sub.index <- sample(1:nrow(train),0.65*nrow(train))
sub.train <- train[sub.index,] #sub train: 25339 obs
validation <- train[-sub.index,] #validation: 13645 obs
Si verifica infine che le proporzioni delle due classi della variabile target siano le stesse su tutti i set di dati creati e che corrispondano a quelle del dataset di partenza (63% non fumatore - 37% fumatore).
round(prop.table(table(sub.train$smoking)),2)
##
## 0 1
## 0.63 0.37
round(prop.table(table(validation$smoking)),2)
##
## 0 1
## 0.63 0.37
round(prop.table(table(test$smoking)),2)
##
## 0 1
## 0.63 0.37
Tutto è correttamente proporzionato, con un leggero sbilanciamento delle classi. Si procede quindi con la fase di pre-processing.
La seguente analisi esplorativa prenderà in considerazione solo il sub.train. Per prima cosa si ispeziona la presenza di valori mancanti.
sum(is.na(sub.train)) #sono presenti missing values?
## [1] 0
md.pattern(sub.train, rotate.names=TRUE)
## /\ /\
## { `---' }
## { O O }
## ==> V <== No need for mice. This data set is completely observed.
## \ \|/ /
## `-----'
## gender age height.cm. weight.kg. waist.cm. hearing.left. hearing.right.
## 25339 1 1 1 1 1 1 1
## 0 0 0 0 0 0 0
## systolic relaxation fasting.blood.sugar Cholesterol triglyceride HDL LDL
## 25339 1 1 1 1 1 1 1
## 0 0 0 0 0 0 0
## hemoglobin serum.creatinine AST ALT Gtp dental.caries tartar smoking
## 25339 1 1 1 1 1 1 1 1 0
## 0 0 0 0 0 0 0 0 0
sub.train.numeric <- sub.train%>% dplyr::select(-c(gender, hearing.left., hearing.right., dental.caries, tartar, smoking))
summary(sub.train.numeric) #statistiche principali per le variabili quantitative
## age height.cm. weight.kg. waist.cm.
## Min. :20.0 Min. :135.0 Min. : 30.00 Min. : 51.00
## 1st Qu.:40.0 1st Qu.:160.0 1st Qu.: 55.00 1st Qu.: 76.00
## Median :40.0 Median :165.0 Median : 65.00 Median : 82.00
## Mean :44.2 Mean :164.7 Mean : 65.86 Mean : 82.03
## 3rd Qu.:55.0 3rd Qu.:170.0 3rd Qu.: 75.00 3rd Qu.: 88.00
## Max. :85.0 Max. :190.0 Max. :130.00 Max. :128.00
## systolic relaxation fasting.blood.sugar Cholesterol
## Min. : 80.0 Min. : 40.00 Min. : 46.00 Min. : 55.0
## 1st Qu.:112.0 1st Qu.: 70.00 1st Qu.: 89.00 1st Qu.:172.0
## Median :120.0 Median : 76.00 Median : 96.00 Median :195.0
## Mean :121.5 Mean : 75.98 Mean : 99.29 Mean :196.8
## 3rd Qu.:130.0 3rd Qu.: 82.00 3rd Qu.:104.00 3rd Qu.:219.0
## Max. :213.0 Max. :146.00 Max. :505.00 Max. :441.0
## triglyceride HDL LDL hemoglobin
## Min. : 11.0 Min. : 4.00 Min. : 1 Min. : 4.90
## 1st Qu.: 74.0 1st Qu.: 47.00 1st Qu.: 92 1st Qu.:13.60
## Median :107.0 Median : 55.00 Median : 113 Median :14.80
## Mean :126.3 Mean : 57.22 Mean : 115 Mean :14.62
## 3rd Qu.:159.0 3rd Qu.: 65.00 3rd Qu.: 136 3rd Qu.:15.70
## Max. :548.0 Max. :618.00 Max. :1660 Max. :20.90
## serum.creatinine AST ALT Gtp
## Min. : 0.1000 Min. : 6.00 Min. : 2.00 Min. : 3.00
## 1st Qu.: 0.8000 1st Qu.: 19.00 1st Qu.: 15.00 1st Qu.: 17.00
## Median : 0.9000 Median : 23.00 Median : 21.00 Median : 25.00
## Mean : 0.8839 Mean : 26.28 Mean : 27.17 Mean : 39.93
## 3rd Qu.: 1.0000 3rd Qu.: 28.00 3rd Qu.: 31.00 3rd Qu.: 43.00
## Max. :10.0000 Max. :1311.00 Max. :2062.00 Max. :999.00
Il dataset non contiene missing. E’ necessario controllare anche la presenza di eventuali outliers. Per fare ciò si riporta una tabella contenente massimo e minimo delle variabili quantitative e il range di variazione ottimale di esse (valori per una persona sana).
Train:
| Variabile | Minimo | Massimo | Range di variazione ottimale |
|---|---|---|---|
| Systolic | 80 | 213 | 110-140 mm/hg |
| Relaxation | 40 | 146 | 70-90 mm/hg |
| Fasting blood sugar | 46 | 505 | 70-140 mg/dl |
| Cholesterol | 55 | 441 | <200 mg/dl |
| Triglyceride | 11 | 548 | <200 mg/dl |
| HDL | 4 | 618 | >40 mg/dl |
| LDL | 1 | 1660 | <130 mg/dl |
| Hemoglobin | 4.9 | 20.9 | 12-18 mg/dl |
| Serum creatinine | 0.1 | 10 | 0.8-1.2 mg/dl |
| AST | 6 | 1311 | 8-60 U/l |
| ALT | 1 | 2062 | 7-55 U/l |
| γ-GTP | 3 | 999 | 2-50 U/l |
Essendo misurazioni mediche è normale che molte osservazioni siano fuori dal range di valori ottimale: ad esempio, un valore di 455 per il colesterolo è molto elevato e pericoloso per la salute, ma non per questo è da considerarsi un valore non plausibile. La presenza di outlier è quindi intrinseca nel tipo di variabili rilevate, inoltre la numerosità dei valori più estremi rispetto al totale delle osservazioni è irrisoria: queste considerazioni hanno portato alla conclusione di non rimuoverli, non causando particolari distorsioni sulle analisi.
correlazione <- cor(sub.train.numeric)
heatmaply::heatmaply_cor(correlazione,
cellnote = round(correlazione,2), cellnote_textposition = "middle center",
cellnote_size = 8,
show_dendrogram = c(FALSE, FALSE))
Si notano alcune correlazioni significative:
Le prime due coppie hanno una correlazione molto elevata, il che giustifica un’ulteriore selezione di variabili. In particolare, si mantengono AST (enzima più generico) e weight (waist fornisce un’informazione ridondante).
train <- train%>% dplyr::select(-c(waist.cm., ALT))
sub.train <- sub.train%>% dplyr::select(-c(waist.cm., ALT))
sub.train.numeric <- sub.train.numeric%>% dplyr::select(-c(waist.cm., ALT))
validation <- validation%>% dplyr::select(-c(waist.cm., ALT))
test <- test%>% dplyr::select(-c(waist.cm., ALT))
Si procede ora con una rappresentazione grafica delle esplicative attraverso dei boxplot condizionati alle classi della variabile target.
boxplot <-list()
variabili<- colnames(sub.train.numeric)
for (i in variabili){
boxplot[[i]] <- ggplot(sub.train, aes_string(x = "smoking", y = i, col = "smoking", fill = "smoking")) +
geom_boxplot(alpha = 0.2) +
theme(legend.position = "none") +
scale_color_manual(values = c("blue", "red"))
scale_fill_manual(values = c("blue", "red"))
}
library(gridExtra)
grid.arrange(grobs=boxplot, nrow=4, ncol=4)
In alcuni boxplot si notano alcuni outlier: in HDL, ad esempio, c’è solo un valore completamente fuori dal range. Si tratta comunque di un osservazione su un sub.train con una numerosità di 25000, per questo non influenzerà le analisi e non verrà quindi eliminato. Si nota inoltre come alcune variabili abbiano una distribuzione fortemente asimmetrica, evidenziata anche dagli istogrammi di esse: per questo è opportuno effettuare una trasformazione logaritmica.
par(mfrow=c(1,2))
hist(sub.train$fasting.blood.sugar, main = 'fasting.blood.sugar')
sub.train$fasting.blood.sugar<-log(sub.train$fasting.blood.sugar)
hist(sub.train$fasting.blood.sugar, main = 'log(fasting.blood.sugar)')
hist(sub.train$triglyceride, main = 'triglyceride')
sub.train$triglyceride<-log(sub.train$triglyceride)
hist(sub.train$triglyceride, main = 'log(triglyceride)')
hist(sub.train$HDL, main = 'HDL')
sub.train$HDL<-log(sub.train$HDL)
hist(sub.train$HDL, main = 'log(HDL)')
hist(sub.train$LDL, main = 'LDL')
sub.train$LDL<-log(sub.train$LDL)
hist(sub.train$LDL, main = 'log(LDL)')
hist(sub.train$serum.creatinine, main = 'serum.creatinine')
sub.train$serum.creatinine<-log(sub.train$serum.creatinine)
hist(sub.train$serum.creatinine, main = 'log(serum.creatinine)')
hist(sub.train$AST, main = 'AST')
sub.train$AST<-log(sub.train$AST)
hist(sub.train$AST, main = 'log(AST)')
hist(sub.train$Gtp, main = 'Gtp')
sub.train$Gtp<-log(sub.train$Gtp)
hist(sub.train$Gtp, main = 'log(Gtp)')
par(mfrow=c(1,1))
Confrontando gli istogrammi delle variabili originali e delle trasformate logaritmiche si nota come la distribuzione sia ora più simmetrica.
La normalizzazione è altamente influenzata da quelli che potrebbero essere dei valori anomali nel sub.train, andando a genererare distorsioni nel validation e successivamente nel test, per questo motivo si è scelto di non utilizzarla. La scelta è ricaduta tra due tecniche, la prima è una standardizzazione e la seconda una tecnica costruita ad hoc andando a sottrarre al valore il primo quartile e dividendo il tutto per lo scarto interquartile. Si è preferito usare la prima: nonostante la media sia influenzata dagli outlier, la numerosità elevata mitiga questo effetto e si ottiene una variabile con distribuzione nota.
#Creazione matrice per tenerne traccia di min e max di ogni variabile
matrix_indicators <- matrix(0, nrow=(dim(sub.train.numeric)[2]), ncol=2)
colnames(matrix_indicators) <- c("media", "dv")
sub.train.numeric$fasting.blood.sugar<-(sub.train$fasting.blood.sugar)
sub.train.numeric$triglyceride<-(sub.train$triglyceride)
sub.train.numeric$HDL<-(sub.train$HDL)
sub.train.numeric$LDL<-(sub.train$LDL)
sub.train.numeric$serum.creatinine<-(sub.train$serum.creatinine)
sub.train.numeric$AST<-(sub.train$AST)
sub.train.numeric$Gtp<-(sub.train$Gtp)
#Calcolo media e deviazione standard
for (i in 1:nrow(matrix_indicators)){
matrix_indicators[i, "media"] <- mean(sub.train.numeric[,i])
matrix_indicators[i, "dv"] <- sd(sub.train.numeric[,i])
}
#Standardizzazione
for (i in 1:nrow(matrix_indicators)){
sub.train.numeric[, i] <- (sub.train.numeric[, i] - matrix_indicators[i, "media"])/(matrix_indicators[i, "dv"]) }
#Aggiornamento sub.train completo (con anche le qualitative)
sub.train <- cbind(sub.train.numeric, sub.train$gender, sub.train$hearing.left., sub.train$hearing.right.,
sub.train$dental.caries, sub.train$tartar, sub.train$smoking)
colnames(sub.train)[15:20] <- c('gender','hearing.left.','hearing.right.','dental.caries','tartar','smoking')
Viene svolta la analisi per componenti principali al fine sia di individuare le variabili che più influenzano i pesi delle componenti, sia per graficare la distribuzione della classificazione con le due variabili più legate alle prime due componenti principali. Queste sono trigliceridi e Gtp: verranno tenute in considerazione solo a fini grafici, mentre i modelli verranno stimati con il set di variabili completo.
pca <- princomp(sub.train.numeric[,-c(1:3)])
summary(pca)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 1.7161246 1.3626498 1.1958172 1.0531326 1.03452262
## Proportion of Variance 0.2677454 0.1688080 0.1300032 0.1008302 0.09729812
## Cumulative Proportion 0.2677454 0.4365534 0.5665566 0.6673868 0.76468490
## Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
## Standard deviation 0.90435635 0.76818257 0.72088586 0.61261831 0.4849610
## Proportion of Variance 0.07435388 0.05364798 0.04724518 0.03411964 0.0213815
## Cumulative Proportion 0.83903878 0.89268676 0.93993194 0.97405157 0.9954331
## Comp.11
## Standard deviation 0.224129940
## Proportion of Variance 0.004566928
## Cumulative Proportion 1.000000000
#le prime due componenti
sort(pca$loadings[,1]) #Gtp
## HDL LDL Cholesterol fasting.blood.sugar
## -0.24483128 0.09973841 0.15332011 0.23547218
## serum.creatinine AST systolic relaxation
## 0.24401486 0.28763844 0.34607437 0.35833805
## hemoglobin triglyceride Gtp
## 0.35970767 0.38665206 0.42639597
sort(pca$loadings[,2]) #Cholesterol
## Cholesterol LDL HDL triglyceride
## -0.69058207 -0.67570481 -0.18595675 -0.01900431
## AST Gtp relaxation hemoglobin
## 0.02604915 0.04414287 0.04468465 0.06289490
## systolic serum.creatinine fasting.blood.sugar
## 0.07491514 0.08022087 0.10464722
ggplot(sub.train, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point(alpha=0.2)+
stat_ellipse(type="norm")+
scale_colour_manual(values=c("darkgreen", "red"))+
theme_bw()+
labs(title="Classificazione originale")+
xlab("trasformazione Gtp")+
ylab("trasformazione Colesterolo")
Ogni modello verrà costruito sul sub train e testato nel validation, verrà creata una matrice che salva le seguenti metriche: F1 score, accuracy e recall(sensitivity). L’F1 score risulta particolarmente vantaggiosa in questo caso per due motivi. In primis è utile quando ci sono costi diversi di falsi positivi e falsi negativi: in questo caso, per un’assicurazione è peggio predire come non fumatori persone che in realtà lo sono, in quanto i costi calcolati sarebbero minore di quelli potenzialmente necessari. Inoltre le classi non sono perfettamente bilanciate (63%-37%), quindi limitarsi all’accuracy potrebbe dare risultati imprecisi.
matrice_ris<-matrix(ncol=4)
colnames(matrice_ris)<-c("modello", "F1", "accuracy", "recall")
L’intera analisi dei dati, data anche la numerosità del dataset, verrà effettuata valutando le prestazioni dei modelli costruiti sul sub.train nel validation. Si confronteranno successivamente i risultati ottenuti nel validation al fine di scegliere un unico modello da testare sul test set. L’obiettivo è simulare una situazione reale in cui il test set rappresenta i soggetti su cui usare il modello per predire la loro condizione di fumatori.
Il primo modello testato è un albero di classificazione. I limiti di questo modello sono in primis a livello computazionale. Difatti vista la grande dimensionalità del dataset il numero di nodi che il modello va a creare è estremamente esiguo (3). Si è provato a rieseguire l’analisi sia riducendo la dimensionalità che usando le variabili senza le trasformazioni ma non sono stati ottenuti risultati migliori. Inoltre non è stata necessaria l’applicazione della potatura dell’albero in quanto già di dimensioni ridotte.
albero <- tree(smoking~., data=sub.train)
summary(albero)
##
## Classification tree:
## tree(formula = smoking ~ ., data = sub.train)
## Variables actually used in tree construction:
## [1] "gender" "Gtp"
## Number of terminal nodes: 3
## Residual mean deviance: 0.9717 = 24620 / 25340
## Misclassification error rate: 0.2891 = 7325 / 25339
plot(albero)
text(albero, pretty = 0)
L’albero ottenuto classifica inizialmente sul gender: se il genere è femminile la classificazione prevista è non fumatrici, se il genere è maschile invece viene effettuata una discriminazione in base al gtp. Come è immaginabile un modello così semplice risulta poco informativo, tanto che la misclassificazione è quasi del 29%.
albero
## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 25339 33290 0 ( 0.63396 0.36604 )
## 2) gender: F 9259 3063 0 ( 0.96079 0.03921 ) *
## 3) gender: M 16080 22100 1 ( 0.44577 0.55423 )
## 6) Gtp < 0.606893 10756 14910 0 ( 0.50958 0.49042 ) *
## 7) Gtp > 0.606893 5324 6650 1 ( 0.31687 0.68313 ) *
Nello split basato sul Gtp la proporzione delle classi è quasi identica (55%-45%), questo rende la discriminazione legata al training set. Inoltre è interessante notare la proporzione di donne fumatrici (meno del 4%).
Prima di testare il modello sul validation, è opportuno e necessario effettuare su esso le trasformazioni svolte nel sub train.
summary(validation[,-20]) #no missing values
## gender age height.cm. weight.kg. hearing.left.
## F:4960 Min. :20.00 Min. :135.0 Min. : 30.00 1:13292
## M:8685 1st Qu.:40.00 1st Qu.:160.0 1st Qu.: 55.00 2: 353
## Median :40.00 Median :165.0 Median : 65.00
## Mean :44.14 Mean :164.6 Mean : 65.82
## 3rd Qu.:55.00 3rd Qu.:170.0 3rd Qu.: 75.00
## Max. :85.00 Max. :190.0 Max. :135.00
## hearing.right. systolic relaxation fasting.blood.sugar
## 1:13286 Min. : 79.0 Min. : 46.00 Min. : 46.0
## 2: 359 1st Qu.:112.0 1st Qu.: 70.00 1st Qu.: 89.0
## Median :120.0 Median : 76.00 Median : 96.0
## Mean :121.6 Mean : 76.12 Mean : 99.5
## 3rd Qu.:130.0 3rd Qu.: 82.00 3rd Qu.:104.0
## Max. :240.0 Max. :146.00 Max. :398.0
## Cholesterol triglyceride HDL LDL
## Min. : 77.0 Min. : 16.0 Min. : 18.00 Min. : 1.0
## 1st Qu.:172.0 1st Qu.: 74.0 1st Qu.: 47.00 1st Qu.: 91.0
## Median :195.0 Median :108.0 Median : 56.00 Median : 112.0
## Mean :196.7 Mean :126.6 Mean : 57.41 Mean : 114.9
## 3rd Qu.:220.0 3rd Qu.:161.0 3rd Qu.: 66.00 3rd Qu.: 136.0
## Max. :445.0 Max. :999.0 Max. :155.00 Max. :1860.0
## hemoglobin serum.creatinine AST Gtp
## Min. : 4.90 Min. : 0.1000 Min. : 6.00 Min. : 1.00
## 1st Qu.:13.60 1st Qu.: 0.8000 1st Qu.: 19.00 1st Qu.: 17.00
## Median :14.80 Median : 0.9000 Median : 23.00 Median : 26.00
## Mean :14.62 Mean : 0.8854 Mean : 26.23 Mean : 40.21
## 3rd Qu.:15.70 3rd Qu.: 1.0000 3rd Qu.: 29.00 3rd Qu.: 43.00
## Max. :21.10 Max. :11.6000 Max. :981.00 Max. :976.00
## dental.caries tartar
## 0:10736 N:5986
## 1: 2909 Y:7659
##
##
##
##
#Applicazione trasformate logaritmiche
validation$fasting.blood.sugar<-log(validation$fasting.blood.sugar)
validation$triglyceride<-log(validation$triglyceride)
validation$HDL<-log(validation$HDL)
validation$LDL<-log(validation$LDL)
validation$serum.creatinine<-log(validation$serum.creatinine)
validation$AST<-log(validation$AST)
validation$Gtp<-log(validation$Gtp)
#Standardizzazione
validation.numeric <- validation[,-c(1,5,6,18,19,20)]
for (i in 1:nrow(matrix_indicators)){
validation.numeric[, i] <- (validation.numeric[, i] - matrix_indicators[i, "media"])/(matrix_indicators[i, "dv"])
}
validation <- cbind(validation.numeric, validation$gender, validation$hearing.left., validation$hearing.right.,
validation$dental.caries, validation$tartar, validation$smoking)
colnames(validation)[15:20] <- c('gender','hearing.left.','hearing.right.','dental.caries','tartar','smoking')
Si valuta l’efficienza sul validation.
tree.pred <- predict(albero , validation[,-20], type = "class")
cfalbero<-confusionMatrix(tree.pred , validation$smoking, positive = "1", mode="everything")
## Registered S3 methods overwritten by 'proxy':
## method from
## print.registry_field registry
## print.registry_entry registry
cfalbero
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 7731 3050
## 1 909 1955
##
## Accuracy : 0.7099
## 95% CI : (0.7022, 0.7175)
## No Information Rate : 0.6332
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.3136
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.3906
## Specificity : 0.8948
## Pos Pred Value : 0.6826
## Neg Pred Value : 0.7171
## Precision : 0.6826
## Recall : 0.3906
## F1 : 0.4969
## Prevalence : 0.3668
## Detection Rate : 0.1433
## Detection Prevalence : 0.2099
## Balanced Accuracy : 0.6427
##
## 'Positive' Class : 1
##
L’accuracy del 71% è fuorviante, infatti la classificazione dei non fumatori è efficiente (specificità quasi del 90%), ma facendo riferimento ai fumatori il modello risulta altamente inefficiente (sensibilità al di sotto del 40%).
#Salva i risultati ottenuti nella matrice
matrice_ris[1,]<-c("albero decionale", cfalbero$byClass["F1"], cfalbero$overall["Accuracy"],cfalbero$byClass["Sensitivity"])
Ora si provano a graficare i risultati ottenuti con la classificazione e le vere label, evidenziando in blu i punti misclassificati.
library(mclust)
## Warning: package 'mclust' was built under R version 4.2.3
## Package 'mclust' version 6.0.1
## Type 'citation("mclust")' for citing this R package in publications.
##
## Attaching package: 'mclust'
## The following object is masked from 'package:purrr':
##
## map
misclassificatealbero<-(classError(tree.pred, validation$smoking)$misclassified)
plotv<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"))+
theme_bw()+
labs("Classificazione originale vere label")+
theme_bw()
plotc<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
geom_point(data=validation[misclassificatealbero,], col="blue", alpha=0.5)+
theme_bw()+
labs("Classificazione secondo il metodo tree")+
theme_bw()
grid.arrange(plotv, plotc, ncol = 2,
top = "Etichette vere vs classificazione con un albero")
Il blu indica una classificazione errata. Nel grafico si nota come l’albero abbia una grossa difficoltà a classificare principalmente le unità “centrali”.
Il modello più restrittivo in termini di ipotesi sulle covariate è quello basato su LDA/QDA. Questi metodi impongono in primis che sia verificata la distribuzione normale delle variabili. Per avere una visione immediata vengono qui presentati i QQplot condizionati alla classe; ad essi viene sovraimposta la linea dei quantili ‘teorici’ in modo da avere un’indicazione dello scostamento di quelli della variabile in esame.
Da questa analisi sono escluse le variabili dicotomiche, age, height e weight. Osservando i valori di queste ultime 3 si nota infatti come assumano valori ad intervalli di 5 unità: le rilevazioni sono quindi state arrotondate. Si potrebbe infatti pensare di ricodificarle utilizzando degli intervalli e renderle variabili factor. Al di fuori di esse rimane comunque un ampio set di variabili, quindi si procede con esse.
#QQplot condizionati alla classe smoking = 1 (fumatori)
par(mfrow = c(2,3))
for(i in variabili[4:9]) {
qqnorm(sub.train[sub.train$smoking == 1, i], main = i)
qqline(sub.train[sub.train$smoking == 1, i], col = 2)
}
par(mfrow = c(2,3))
for(i in variabili[10:14]) {
qqnorm(sub.train[sub.train$smoking == 1, i], main = i)
qqline(sub.train[sub.train$smoking == 1, i], col = 2)
}
#QQplot condizionati alla classe smoking = 0 (non fumatori)
par(mfrow = c(2,3))
for(i in variabili[4:9]) {
qqnorm(sub.train[sub.train$smoking == 0, i], main = i)
qqline(sub.train[sub.train$smoking == 0, i], col = 2)
}
par(mfrow = c(2,3))
for(i in variabili[10:14]) {
qqnorm(sub.train[sub.train$smoking == 0, i], main = i)
qqline(sub.train[sub.train$smoking == 0, i], col = 2)
}
Già dai QQplot si nota come nessuna delle variabili abbia un andamento normale per entrambe le classi. Per avere un’ulteriore conferma si costruiscano le curve di densità condizionate.
plot_density <- list()
for(i in variabili[4:14]){
plot_density[[i]] <- ggplot(sub.train, aes_string(x = i, y = "..density..", col = "smoking")) +
geom_density(aes(y = ..density..)) +
scale_color_manual(values = c("blue", "red")) +
theme(legend.position = "none")
}
do.call(grid.arrange, c(plot_density, nrow = 4))
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
I problemi sono soprattutto sulle code: questo è causato dai numerosi valori estremi che caratterizzano le variabili. Un’ulteriore conferma sarebbe fornita dal test per la normalità di Shapiro-Wilk, non applicabile in quest’analisi a causa della numerosità delle osservazioni molto maggiore rispetto a quella richiesta per lo svolgimento del test (dal punto di vista computazionale).
Dall’analisi fatta a livello qualitativo sui grafici delle densità, alcune variabili come colesterolo o HDL hanno una distribuzione normale. Il problema risiede nella perfetta sovrapposizione delle due curve di densità condizionate, che impediscono ai modelli di identificare in maniera efficace le due classi. E’ stata comunque svolta un’analisi per un’ulteriore conferma, la quale ha dato pessimi risultati. Le uniche variabili che sembrano essere distinte sono i trigliceridi, emoglobina e Gtp ma non presentano una distribuzione normale, potrebbero essere misture di misture. Si sceglie di procedere con metodi ritenuti più opportuni.
Si applica ora la regressione logistica. Prima si stima un modello con tutte le esplicative e successivamente, per cercare di ridurre la dimensionalità, si usa una selezione stepwise basata sul criterio AIC (the lower the better) in modo da rimuovere le variabili non significative.
logit <- glm(smoking ~., data = sub.train, family = binomial)
summary(logit) #AIC: 23405
##
## Call:
## glm(formula = smoking ~ ., family = binomial, data = sub.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5620 -0.8798 -0.2321 0.8924 3.0405
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.213295 0.067805 -47.391 < 2e-16 ***
## age -0.025855 0.018746 -1.379 0.16782
## height.cm. 0.211427 0.029391 7.194 6.31e-13 ***
## weight.kg. -0.220915 0.024554 -8.997 < 2e-16 ***
## systolic -0.197098 0.025913 -7.606 2.83e-14 ***
## relaxation 0.045124 0.025163 1.793 0.07293 .
## fasting.blood.sugar 0.026668 0.017003 1.568 0.11678
## Cholesterol -0.001838 0.050627 -0.036 0.97104
## triglyceride 0.276122 0.028689 9.625 < 2e-16 ***
## HDL -0.026433 0.027175 -0.973 0.33071
## LDL -0.119755 0.045465 -2.634 0.00844 **
## hemoglobin 0.171546 0.025415 6.750 1.48e-11 ***
## serum.creatinine -0.200441 0.021282 -9.418 < 2e-16 ***
## AST -0.226599 0.018803 -12.051 < 2e-16 ***
## Gtp 0.537373 0.022545 23.835 < 2e-16 ***
## genderM 2.992413 0.079223 37.772 < 2e-16 ***
## hearing.left.2 -0.115485 0.120547 -0.958 0.33806
## hearing.right.2 0.092579 0.119591 0.774 0.43886
## dental.caries1 0.310428 0.038648 8.032 9.57e-16 ***
## tartarY 0.306661 0.032797 9.350 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 33286 on 25338 degrees of freedom
## Residual deviance: 23365 on 25319 degrees of freedom
## AIC: 23405
##
## Number of Fisher Scoring iterations: 6
step.logit <- stepAIC(logit, direction = "both", trace = FALSE)
summary(step.logit) #AIC: 23400
##
## Call:
## glm(formula = smoking ~ height.cm. + weight.kg. + systolic +
## relaxation + triglyceride + LDL + hemoglobin + serum.creatinine +
## AST + Gtp + gender + dental.caries + tartar, family = binomial,
## data = sub.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5896 -0.8816 -0.2323 0.8921 3.0394
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.22105 0.06733 -47.838 < 2e-16 ***
## height.cm. 0.21405 0.02795 7.659 1.88e-14 ***
## weight.kg. -0.20590 0.02353 -8.752 < 2e-16 ***
## systolic -0.20045 0.02564 -7.818 5.36e-15 ***
## relaxation 0.04526 0.02512 1.802 0.0716 .
## triglyceride 0.28822 0.01847 15.601 < 2e-16 ***
## LDL -0.12326 0.01634 -7.543 4.60e-14 ***
## hemoglobin 0.17642 0.02494 7.075 1.50e-12 ***
## serum.creatinine -0.20238 0.02122 -9.537 < 2e-16 ***
## AST -0.22901 0.01871 -12.238 < 2e-16 ***
## Gtp 0.53476 0.02208 24.218 < 2e-16 ***
## genderM 3.00038 0.07799 38.473 < 2e-16 ***
## dental.caries1 0.31451 0.03846 8.179 2.87e-16 ***
## tartarY 0.30991 0.03275 9.463 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 33286 on 25338 degrees of freedom
## Residual deviance: 23372 on 25325 degrees of freedom
## AIC: 23400
##
## Number of Fisher Scoring iterations: 6
Le variabili significative mantenute dalla selezione stepwise sono comunque molto numerose. Si procede analizzando la presenza di punti influenti, eliminandoli e ripetendo la stima del modello per vedere se ci sono dei miglioramenti.
influencePlot(step.logit)
## StudRes Hat CookD
## 20340 -1.675502 8.783723e-03 0.0019259620
## 4225 -2.468969 2.086182e-03 0.0029393022
## 11570 1.682581 8.739313e-03 0.0019461947
## 1577 3.040595 6.957062e-05 0.0004990126
## 17822 -2.598540 1.683238e-03 0.0033280060
## 16607 3.004504 8.576728e-05 0.0005507912
sub.train2 <- sub.train[-c(20340,4225,11570,1577,17822,16607),]
logit2 <- glm(smoking ~., data = sub.train2, family = binomial)
summary(logit2) #AIC: 23402
##
## Call:
## glm(formula = smoking ~ ., family = binomial, data = sub.train2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5612 -0.8802 -0.2320 0.8925 3.0398
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.212615 0.067802 -47.382 < 2e-16 ***
## age -0.025509 0.018747 -1.361 0.17362
## height.cm. 0.211562 0.029392 7.198 6.11e-13 ***
## weight.kg. -0.220420 0.024556 -8.976 < 2e-16 ***
## systolic -0.197443 0.025917 -7.618 2.57e-14 ***
## relaxation 0.045226 0.025164 1.797 0.07230 .
## fasting.blood.sugar 0.026474 0.017007 1.557 0.11955
## Cholesterol 0.001999 0.050975 0.039 0.96871
## triglyceride 0.274407 0.028788 9.532 < 2e-16 ***
## HDL -0.028021 0.027269 -1.028 0.30415
## LDL -0.122991 0.045757 -2.688 0.00719 **
## hemoglobin 0.171089 0.025417 6.731 1.68e-11 ***
## serum.creatinine -0.200285 0.021282 -9.411 < 2e-16 ***
## AST -0.226378 0.018804 -12.039 < 2e-16 ***
## Gtp 0.537272 0.022544 23.832 < 2e-16 ***
## genderM 2.992071 0.079222 37.768 < 2e-16 ***
## hearing.left.2 -0.115652 0.120543 -0.959 0.33735
## hearing.right.2 0.092334 0.119589 0.772 0.44005
## dental.caries1 0.310326 0.038650 8.029 9.82e-16 ***
## tartarY 0.306032 0.032799 9.331 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 33277 on 25332 degrees of freedom
## Residual deviance: 23362 on 25313 degrees of freedom
## AIC: 23402
##
## Number of Fisher Scoring iterations: 6
step.logit2 <- stepAIC(logit2, direction = "both", trace = FALSE)
summary(step.logit2) #AIC: 23396
##
## Call:
## glm(formula = smoking ~ height.cm. + weight.kg. + systolic +
## relaxation + triglyceride + LDL + hemoglobin + serum.creatinine +
## AST + Gtp + gender + dental.caries + tartar, family = binomial,
## data = sub.train2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5886 -0.8817 -0.2324 0.8921 3.0390
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.22038 0.06733 -47.831 < 2e-16 ***
## height.cm. 0.21398 0.02795 7.656 1.92e-14 ***
## weight.kg. -0.20542 0.02353 -8.730 < 2e-16 ***
## systolic -0.20081 0.02564 -7.831 4.83e-15 ***
## relaxation 0.04541 0.02512 1.808 0.0706 .
## triglyceride 0.28811 0.01847 15.596 < 2e-16 ***
## LDL -0.12326 0.01634 -7.543 4.60e-14 ***
## hemoglobin 0.17586 0.02494 7.052 1.77e-12 ***
## serum.creatinine -0.20214 0.02122 -9.525 < 2e-16 ***
## AST -0.22868 0.01871 -12.220 < 2e-16 ***
## Gtp 0.53467 0.02208 24.215 < 2e-16 ***
## genderM 3.00001 0.07798 38.470 < 2e-16 ***
## dental.caries1 0.31436 0.03846 8.174 2.98e-16 ***
## tartarY 0.30925 0.03275 9.442 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 33277 on 25332 degrees of freedom
## Residual deviance: 23368 on 25319 degrees of freedom
## AIC: 23396
##
## Number of Fisher Scoring iterations: 6
Il modello risultante comprende le seguenti variabili: height, weight, systolic, relaxation, triglyceride, LDL, hemoglobin, serum.creatinine, AST, Gtp, gender, dental.caries, tartar.
La regressione logistica assume che ci sia una relazione lineare tra le variabili esplicative e il logaritmo dell’odds ratio: è quindi bene testare ciò attraverso un’analisi grafica.
#Calcolo fitted values (focus solo sulle variabili quantitative)
probabilities <- predict(step.logit2, type = "response") #posterior probabilities
variabili.new <- variabili[-(6:7)]
#Costruzione degli odds ratio
supp <- sub.train2[,variabili.new]
supp <- supp %>%
mutate(logit = log(probabilities/(1-probabilities))) %>%
gather(key = "variabili.new", value = "predictor.value", -logit)
#Costruzione dei grafici
ggplot(supp, aes(logit, predictor.value))+
geom_point(size = 0.5, alpha = 0.5) +
geom_smooth(method = "loess") +
theme_bw() +
facet_wrap(~variabili[-(6:7)], scales = "free_y")
## `geom_smooth()` using formula = 'y ~ x'
Come si nota la linearità è in parte rispettata ma l’alta variabilità dei punti non permette un’investigazione più precisa. Inoltre c’è un’ulteriore evidenza del fatto che age, height e weight andrebbero trasformate in dummy individuando più intervalli a cui assegnare le osservazioni.
predict.logit <- predict(step.logit2, validation[, -20], type = "response") #posterior probabilities
predict.logit.class <- ifelse(predict.logit > 0.5, 1, 0) #suddivide le post prob in due classi (0-1) con la soglia decisionale pari a 0.5
conf <- confusionMatrix(factor(validation[, 20]), factor(as.vector(predict.logit.class)), positive = '1', mode='everything')
conf
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 6745 1895
## 1 1505 3500
##
## Accuracy : 0.7508
## 95% CI : (0.7435, 0.7581)
## No Information Rate : 0.6046
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4722
##
## Mcnemar's Test P-Value : 2.536e-11
##
## Sensitivity : 0.6487
## Specificity : 0.8176
## Pos Pred Value : 0.6993
## Neg Pred Value : 0.7807
## Precision : 0.6993
## Recall : 0.6487
## F1 : 0.6731
## Prevalence : 0.3954
## Detection Rate : 0.2565
## Detection Prevalence : 0.3668
## Balanced Accuracy : 0.7332
##
## 'Positive' Class : 1
##
La sensitività è meno del 65%, l’F1 score è circa il 67%. Perciò si provi a modificare la soglia per l’assegnazione al fine di migliorare la classificazione dei fumatori.
predict.logit.class2 <- ifelse(predict.logit > 0.3, 1, 0)
conf1 <- confusionMatrix(factor(as.vector(predict.logit.class2)),factor(validation[, 20]),positive='1',
mode='everything')
conf1
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 5080 334
## 1 3560 4671
##
## Accuracy : 0.7146
## 95% CI : (0.707, 0.7222)
## No Information Rate : 0.6332
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.459
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9333
## Specificity : 0.5880
## Pos Pred Value : 0.5675
## Neg Pred Value : 0.9383
## Precision : 0.5675
## Recall : 0.9333
## F1 : 0.7058
## Prevalence : 0.3668
## Detection Rate : 0.3423
## Detection Prevalence : 0.6032
## Balanced Accuracy : 0.7606
##
## 'Positive' Class : 1
##
Il valore della recall è aumentato nettamente, arrivando al 93%. Anche l’F1 score è migliorato. Diminuendo ulteriormente la soglia migliora la recall ma peggiora la precision, quindi si sceglie di mantenere la soglia a 0.3.
matrice_ris<-rbind(matrice_ris,c("regressione logistica",conf1$byClass["F1"],conf1$overall["Accuracy"],conf1$byClass["Sensitivity"]))
Si graficano quindi i risultati ottenuti.
misclassificatereg<-(classError(predict.logit.class2, validation$smoking)$misclassified)
plotv<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"))+
theme_bw()+
labs("Classificazione originale vere label")+
theme_bw()
plotc<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
geom_point(data=validation[misclassificatealbero,], col="blue", alpha=0.5)+
theme_bw()+
labs("Classificazione secondo la regressione logistica")+
theme_bw()
grid.arrange(plotv, plotc, ncol = 2,
top = "Etichette vere vs classificazione con regressione logistica")
Si valuti ora la curva ROC e si calcoli l’AUC.
roc.logit <- prediction(predict.logit, validation[, 20])
performance.logit <- performance(roc.logit,"tpr","fpr")
auc.logit <- performance(roc.logit, measure = "auc")@y.values
auc.logit
## [[1]]
## [1] 0.8351383
plot(performance.logit, colworize = TRUE, main = "Regressione Logistica")
Il modello di regressione logistica non si adatta bene ai dati per il problema già citato della linearità. Inoltre un problema riscontrato è dovuto alla dimensionalità. Sono infatti state effettuate più prove per ridurre il numero di variabili come: analisi delle componenti principali, eliminazione di variabili più specifiche (es. rimuovere HDL e LDL e mantenere Cholesterol). Questi tentativi hanno però condotto a modelli con un AIC maggiore e senza miglioramenti in termini di performance: questo ha portato alla decisione di quasi tutte le variabili, a discapito di un modello molto ricco.
Si sceglie ora di valutare un modello non parametrico basato sulle distanze: per questo motivo verrano usate solo le variabili quantitative già opportunamente trasformate. La funzione KNN chiede in input la risposta e le covariate in due dataset distinti, sia per il train che per il validation.
X_train_picc<-sub.train%>%dplyr::select(-c(smoking,gender,hearing.left., hearing.right., dental.caries, tartar))#SOLO QUANTITATIVE
Y_train_picc<-sub.train$smoking
#validation
X_validation<-validation%>%dplyr::select(-c(smoking,gender,hearing.left., hearing.right., dental.caries, tartar))
Y_validation<-validation$smoking
Si prova ora a trovare il numero di k vicini ottimale andando a confrontare i diversi valori dell’F1 per ciascun k.
k_to_try = 1:100
err_k = rep(x = 0, times = length(k_to_try))
for (i in seq_along(k_to_try)) {
pred = knn(train = X_train_picc,
test = X_validation, cl = Y_train_picc,
k = k_to_try[i])
err_k[i] = confusionMatrix(pred, Y_validation, positive='1')$byClass["F1"]}
L’esecuzione del ciclo risulta esssere onerosa computazionalmente ma i vantaggi operativi sono notevoli. Tramite un plot dei risultati ottenuti con l’esecuzione del ciclo for, si valutano graficamente le performance del modello al variare del parametro di tuning k.
F1_scores<-matrix(nrow=length(k_to_try), ncol=2)
F1_scores[,1]<-1:100
F1_scores[,2]<-err_k
colnames(F1_scores)<-c("k", "F1")
F1gg<-as.data.frame(F1_scores)
ggplot(data=F1gg, aes(x=k, y=F1))+
geom_point(col="blue")+
geom_line(lty="dotted")+
theme_bw()+
geom_text(data=F1gg[91,], aes(label=round(F1,2)), vjust=-0.8)+
geom_hline(yintercept =F1gg[91,2], col="red", alpha=0.4)+
scale_x_continuous(breaks = c(1,seq(5,100, by=5), 91) )+
labs("f1 score in base al k")
Il miglior k è risultato essere quello uguale a 91. Si ristima il modello sul sub train valutandone la performance nel validation.
#miglior k
best_k<-knn(train = X_train_picc,
test = X_validation, cl = Y_train_picc,
k = 91, prob=T)
conf2 <- confusionMatrix(best_k, Y_validation, positive='1', mode='everything')
conf2
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 6864 1773
## 1 1776 3232
##
## Accuracy : 0.7399
## 95% CI : (0.7325, 0.7473)
## No Information Rate : 0.6332
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.4401
##
## Mcnemar's Test P-Value : 0.9732
##
## Sensitivity : 0.6458
## Specificity : 0.7944
## Pos Pred Value : 0.6454
## Neg Pred Value : 0.7947
## Precision : 0.6454
## Recall : 0.6458
## F1 : 0.6456
## Prevalence : 0.3668
## Detection Rate : 0.2369
## Detection Prevalence : 0.3670
## Balanced Accuracy : 0.7201
##
## 'Positive' Class : 1
##
Il modello ha una sensitività e un F1 score non ottimali, rispettivamente intorno al 65%. Questo è dovuto alla struttura dei cluster in cui i punti sono molto sovrapposti tra le due classi. Una soluzione è quella della modifica delle soglie per far si che l’assegnazione avvenga in maniera più coerente rispetto alla reale proporzione delle classi.
prob.knn <- attributes(best_k)$prob
prop_1<-ifelse(best_k == "0", 1-prob.knn, prob.knn)
classi_soglia_knn<-ifelse(prop_1>0.4, 1, 0)
conf3<-confusionMatrix(as.factor(classi_soglia_knn), validation$smoking, positive = "1", mode="everything")
conf3
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 5676 804
## 1 2964 4201
##
## Accuracy : 0.7239
## 95% CI : (0.7163, 0.7313)
## No Information Rate : 0.6332
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.455
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.8394
## Specificity : 0.6569
## Pos Pred Value : 0.5863
## Neg Pred Value : 0.8759
## Precision : 0.5863
## Recall : 0.8394
## F1 : 0.6904
## Prevalence : 0.3668
## Detection Rate : 0.3079
## Detection Prevalence : 0.5251
## Balanced Accuracy : 0.7482
##
## 'Positive' Class : 1
##
La sensibilità, abbassando la soglia a 0.4, è aumentata fino a quasi l’84% e l’F1 score è diventato del 69%.
library(mclust)
misclassificateknn<-(classError(classi_soglia_knn, validation$smoking)$misclassified)
plotv<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"))+
theme_bw()+
labs("Classificazione originale vere label")+
theme_bw()
plotc<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
geom_point(data=validation[misclassificateknn,], col="blue", alpha=0.5)+
theme_bw()+
labs("Classificazione secondo il metodo knn")+
theme_bw()
grid.arrange(plotv, plotc, ncol = 2,
top = "Etichette vere vs classificazione con 91-nn")
matrice_ris<-rbind(matrice_ris,c("91-NN",conf3$byClass["F1"],conf3$overall["Accuracy"],conf3$byClass["Sensitivity"]))
Si valuti ora la curva ROC e si calcoli l’AUC.
prob.knn <- attributes(best_k)$prob
prob.knn <- 2*ifelse(best_k == "0", 1-prob.knn, prob.knn) - 1
roc.knn <- prediction(prob.knn, validation[, 20])
performance.knn <- performance(roc.knn,"tpr","fpr")
auc.knn <- performance(roc.knn, measure = "auc")@y.values
auc.knn
## [[1]]
## [1] 0.8188348
plot(performance.knn, colworize = TRUE, main = "91-NN")
Il metodo Random Forest si basa sull’idea di combinare i risultati di più alberi decisionali.
Nella prima fase dell’algoritmo vengono addestrati più alberi decisionali a partire da campioni casuali estratti con reinserimento dal dataset di partenza (più training set). In questo modo si riduce l’overfitting: ciascun albero è legato ai suoi dati di training, ma la foresta non lo è.
Per quanto riguarda la costruzione di ciascun albero le features da valutare in ogni nodo sono selezionate casualmente: in questo modo gli alberi saranno meno correlati tra di loro e più diversi, riducendo la varianza e aumentando le prestazioni del modello.
Quando si considera una nuova osservazione basta ottenere la previsione della classe da ogni albero decisionale: la previsione finale è quella più frequente guardando all’intera foresta.
Dato che Random Forest non necessita di particolari ipotesi sulle variabili sono state fatte più prove sia usando il dataset originale che quello trasformato (trasformazioni logaritmiche e standardizzazione). I risultati sono pressochè identici, quindi per coerenza procediamo con le variabili standardizzate.
set.seed(123)
rf <- randomForest(smoking ~., data=sub.train,
importance=TRUE, #l'argomento importance = TRUE fa sì che venga valutata l'importanza relativa di ciascun esplicativa
cutoff=c(0.7,0.3)) #regola la soglia decisionale
#valori previsti sul validation
rf.predict1 <- predict(rf, validation[,-20], type='prob')
rf.predict1.class <- predict(rf, validation[,-20])
conf4 <- confusionMatrix(rf.predict1.class, validation$smoking, positive='1', mode='everything')
conf4
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 5660 248
## 1 2980 4757
##
## Accuracy : 0.7634
## 95% CI : (0.7562, 0.7705)
## No Information Rate : 0.6332
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.5432
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9504
## Specificity : 0.6551
## Pos Pred Value : 0.6148
## Neg Pred Value : 0.9580
## Precision : 0.6148
## Recall : 0.9504
## F1 : 0.7467
## Prevalence : 0.3668
## Detection Rate : 0.3486
## Detection Prevalence : 0.5670
## Balanced Accuracy : 0.8028
##
## 'Positive' Class : 1
##
Il valore della sensitivity è ottimo: 95%. I risultati sono molto buoni anche in termini di F1 score. Il fatto di aver spostato la soglia è andato a discapito della specificity che però, per questa analisi, è meno rilevante.
misclassificaterf<-(classError(rf.predict1.class, validation$smoking)$misclassified)
plotv<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"))+
theme_bw()+
labs("Classificazione originale vere label")+
theme_bw()
plotc<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
geom_point(data=validation[misclassificaterf,], col="blue", alpha=0.5)+
theme_bw()+
labs("Classificazione secondo il metodo random forest")+
theme_bw()
grid.arrange(plotv, plotc, ncol = 2,
top = "Etichette vere vs classificazione con random forest")
matrice_ris<-rbind(matrice_ris,c("random forest",conf4$byClass["F1"],conf4$overall["Accuracy"],conf4$byClass["Sensitivity"]))
#curva ROC e AUC
roc.rf <- prediction(rf.predict1[,2], validation[, 20])
performance.rf <- performance(roc.rf,"tpr","fpr")
auc.rf <- performance(roc.rf, measure = "auc")@y.values
auc.rf
## [[1]]
## [1] 0.8882169
plot(performance.rf, colworize = TRUE, main = "Random Forest")
Il coefficiente AUC è di 0.89: si tratta di un ottimo risultato, considerando che un classificatore perfetto ha AUC=1.
Si proceda ora andando a investigare l’importanza relativa di ciascuna delle variabili esplicative:
importance <- importance(rf)
var.imp <- tibble(variabile=row.names(importance),
importanza = importance[,'MeanDecreaseGini'])
ggplot(var.imp, aes(x=reorder(variabile,desc(importanza)), y=importanza, fill=importanza))+
geom_col()+
labs(x="Predittori", y ="Importanza", title="Importanza dei predittori")+
theme(axis.text.x = element_text(angle=90))
L’importanza delle covariate è tanto maggiore quanto è il valore dell’indice di impurità di Gini. Le variabili più importanti per la classificazione sono genere, gtp, emoglobina e in generale le variabili legate al colesterolo. Minor importanza hanno quelle legate all’udito e all’igiene dentale.
Si confrontino innanzitutto i risultati ottenuti dai modelli precedenti in termini di metriche.
matrice_ris<-as.data.frame(matrice_ris)
matrice_ris$F1<-as.numeric(matrice_ris$F1)
matrice_ris$accuracy<-as.numeric(matrice_ris$accuracy)
matrice_ris$recall<-as.numeric(matrice_ris$recall)
matrice_ris
## modello F1 accuracy recall
## 1 albero decionale 0.4968865 0.7098571 0.3906094
## 2 regressione logistica 0.7058024 0.7146207 0.9332667
## 3 91-NN 0.6903862 0.7238549 0.8393606
## 4 random forest 0.7466646 0.7634298 0.9504496
plot.a<-ggplot(matrice_ris, aes(x=modello,y=F1,col=modello))+
geom_point()+
geom_text(aes(label=round(F1,4)),col='black',vjust=1.5,size=3)+
scale_y_continuous(limits=c(0,1))+
theme_bw()+
xlab('')+
theme(legend.position = 'top')
plot.b<-ggplot(matrice_ris, aes(x=modello,y=accuracy,col=modello))+
geom_point()+
geom_text(aes(label=round(accuracy,4)),col='black',vjust=1.5,size=3)+
scale_y_continuous(limits=c(0,1))+
theme_bw()+
xlab('')+
theme(legend.position = 'none')
plot.c<-ggplot(matrice_ris, aes(x=modello,y=recall,col=modello))+
geom_point()+
geom_text(aes(label=round(recall,4)),col='black',vjust=1.5,size=3)+
scale_y_continuous(limits=c(0,1))+
theme_bw()+
theme(legend.position = 'none')
grid.arrange(plot.a,plot.b,plot.c, nrow = 3, top='Confronto metriche sui vari modelli',heights=c(21.5,15,15))
Come si nota dal grafico il modello che performa meglio in termini di tutte le metriche è random forest. Il peggiore invece, come si era già riscontrato, è l’albero decisionale. L’accuracy non è la metrica più adeguata per questo problema, infatti dà uguale peso agli errori di classificazione. In questo caso invece è bene minimizzare il tasso di falsi negativi (non fumatori) e quindi massimizzare la recall: questa metrica è infatti molto più utile agli scopi di questa analisi ed è anche quella che discrimina maggiormente tra i metodi, avvalorando il fatto che random forest sia quello migliore. L’F1 combina recall e precisione, che misura la qualità della classificazione della classe di interesse (fumatori). Anche questa assume il valore più elevato con il metodo random forest.
Il fatto che quest’ultimo sia il classificatore migliore è visibile anche confrontando le curve di ROC.
par(mfrow=c(1,3))
plot(performance.logit, colworize = TRUE, main = "Regressione Logistica")
plot(performance.knn, colworize = TRUE, main = "91-KNN")
plot(performance.rf, colworize = TRUE, main = "Random Forest")
par(mfrow=c(1,1))
Si va quindi a valutare la performance di random forest sul test set. In primis si uniscano sub.train e validation e si trasformi il test set.
summary(test[,-20]) #no missing values
## gender age height.cm. weight.kg. hearing.left.
## F: 6072 Min. :20.00 Min. :130.0 Min. : 35.00 1:16280
## M:10636 1st Qu.:40.00 1st Qu.:160.0 1st Qu.: 55.00 2: 428
## Median :40.00 Median :165.0 Median : 65.00
## Mean :44.19 Mean :164.6 Mean : 65.92
## 3rd Qu.:55.00 3rd Qu.:170.0 3rd Qu.: 75.00
## Max. :85.00 Max. :190.0 Max. :125.00
## hearing.right. systolic relaxation fasting.blood.sugar
## 1:16271 Min. : 71.0 Min. : 40.00 Min. : 48.0
## 2: 437 1st Qu.:112.0 1st Qu.: 70.00 1st Qu.: 89.0
## Median :120.0 Median : 76.00 Median : 96.0
## Mean :121.4 Mean : 75.96 Mean : 99.2
## 3rd Qu.:130.0 3rd Qu.: 82.00 3rd Qu.:104.0
## Max. :220.0 Max. :146.00 Max. :475.0
## Cholesterol triglyceride HDL LDL hemoglobin
## Min. : 77.0 Min. : 8.0 Min. : 4.0 Min. : 1 Min. : 5.00
## 1st Qu.:173.0 1st Qu.: 75.0 1st Qu.: 47.0 1st Qu.: 92 1st Qu.:13.60
## Median :195.0 Median :109.0 Median : 55.0 Median : 113 Median :14.80
## Mean :197.3 Mean :127.2 Mean : 57.3 Mean : 115 Mean :14.63
## 3rd Qu.:220.0 3rd Qu.:160.0 3rd Qu.: 66.0 3rd Qu.: 136 3rd Qu.:15.80
## Max. :393.0 Max. :399.0 Max. :359.0 Max. :1660 Max. :20.90
## serum.creatinine AST Gtp dental.caries tartar
## Min. : 0.1000 Min. : 7 Min. : 3.00 0:13173 N:7457
## 1st Qu.: 0.8000 1st Qu.: 19 1st Qu.: 17.00 1: 3535 Y:9251
## Median : 0.9000 Median : 23 Median : 26.00
## Mean : 0.8889 Mean : 26 Mean : 39.77
## 3rd Qu.: 1.0000 3rd Qu.: 28 3rd Qu.: 44.00
## Max. :10.3000 Max. :778 Max. :933.00
#Applicazione trasformate logaritmiche
test$fasting.blood.sugar<-log(test$fasting.blood.sugar)
test$triglyceride<-log(test$triglyceride)
test$HDL<-log(test$HDL)
test$LDL<-log(test$LDL)
test$serum.creatinine<-log(test$serum.creatinine)
test$AST<-log(test$AST)
test$Gtp<-log(test$Gtp)
#Standardizzazione
test.numeric <- test[,-c(1,5,6,18,19,20)]
for (i in 1:nrow(matrix_indicators)){
test.numeric[, i] <- (test.numeric[, i] - matrix_indicators[i, "media"])/(matrix_indicators[i, "dv"])
}
test <- cbind(test.numeric, test$gender, test$hearing.left., test$hearing.right.,
test$dental.caries, test$tartar, test$smoking)
colnames(test)[15:20] <- c('gender','hearing.left.','hearing.right.','dental.caries','tartar','smoking')
#Unione sub.train e validation
train.finale <- rbind(sub.train,validation)
Si riallena il modello sul train completo e si valuta sul test.
set.seed(123)
rf.finale <- randomForest(smoking ~., data=train.finale, cutoff=c(0.7,0.3))
rf.pred.class.finale <- predict(rf.finale, test[,-20])
conf5 <- confusionMatrix(rf.pred.class.finale, test$smoking, positive='1', mode='everything')
conf5
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 7311 269
## 1 3222 5906
##
## Accuracy : 0.7911
## 95% CI : (0.7848, 0.7972)
## No Information Rate : 0.6304
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.592
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9564
## Specificity : 0.6941
## Pos Pred Value : 0.6470
## Neg Pred Value : 0.9645
## Precision : 0.6470
## Recall : 0.9564
## F1 : 0.7719
## Prevalence : 0.3696
## Detection Rate : 0.3535
## Detection Prevalence : 0.5463
## Balanced Accuracy : 0.8253
##
## 'Positive' Class : 1
##
La sensitivity è superiore al 95% e l’F1 score è del 77%, evidenziando un’ottima capacità previsiva da parte del modello scelto per la classe di riferimento. L’accuracy è comunque vicina all’80%, evidenziando una capacità del modello di classificare correttamente su entrambe le classi.
Si grafichi la classificazione con il modello scelto evidenziando le osservazioni misclassificate.
misclassificatetest<-(classError(rf.pred.class.finale, test$smoking)$misclassified)
plotv<-ggplot(test, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"))+
labs(title="Classificazione originale vere label nel test set")+
theme_bw()
plotv
plotc<-ggplot(test, aes(x=Gtp, y=Cholesterol, col=smoking))+
geom_point()+
scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
geom_point(data=test[misclassificatetest,], col="blue", alpha=0.5)+
labs(title="Classificazione secondo il metodo random forest nel test set")+
theme_bw()
plotc
I problemi della classificazione rimangono nella parte centrale. A causa dell’elevato numero di osservazioni la rappresentazione è poco informativa: le unità misclassificate sono in tutto 3491 su un totale di 16708. Inoltre, al fine del problema in analisi la percentuale di misclassificate nella classe di fumatori è ancora più bassa.
Questa analisi si conclude quindi con la scelta di una modellistica che possa essere utile alla previsione dello stato di fumatore (o non fumatore) in un caso assicurativo, in cui i costi associati ai diversi tipi di errore incidono in modo rilevante.
Fonti: